Relativistic Breit–Wigner (rel_breitwigner) distribution

The relativistic Breit–Wigner distribution is a continuous distribution on \([0,\infty)\) used in high‑energy physics to model resonances (unstable particles) via the uncertainty in their reconstructed invariant mass.

Compared to the (non‑relativistic) Breit–Wigner / Cauchy form, the relativistic version bakes in the dependence on \(x^2\) that arises naturally when modeling resonance behavior in relativistic kinematics.

Learning goals

• Understand what rel_breitwigner models and how it relates to cauchy / Breit–Wigner.

• Work with the PDF and CDF (including a numerically robust NumPy implementation).

• Compute mean/variance (and understand why higher moments diverge).

• Sample from scratch (NumPy-only) via numerical inverse CDF.

• Use scipy.stats.rel_breitwigner for evaluation, simulation, and fitting.

import platform

import numpy as np

import plotly.graph_objects as go
import os
import plotly.io as pio
from plotly.subplots import make_subplots

import scipy
from scipy import integrate, optimize
from scipy.stats import cauchy, rel_breitwigner

pio.templates.default = "plotly_white"
pio.renderers.default = os.environ.get("PLOTLY_RENDERER", "notebook")

np.set_printoptions(precision=6, suppress=True)
rng = np.random.default_rng(7)

print("Python", platform.python_version())
print("NumPy", np.__version__)
print("SciPy", scipy.__version__)

Python 3.12.9
NumPy 1.26.2
SciPy 1.15.0

1) Title & Classification

• Name: rel_breitwigner

• Type: continuous distribution

• Support (standardized): \(x \in [0, \infty)\)

• Parameter space:

  • shape: \(\rho > 0\)

  • location: \(\text{loc} \in \mathbb{R}\)

  • scale: \(\text{scale} > 0\)

We write (standardized form):

\[X \sim \mathrm{RelBW}(\rho).\]

SciPy uses the name rel_breitwigner and implements the usual location/scale transform:

\[Y = \text{loc} + \text{scale}\,X.\]

2) Intuition & Motivation

2.1 What it models

In collider experiments, unstable particles (resonances) are reconstructed from their decay products. The reconstructed invariant mass is not a single number: it fluctuates due to the particle’s finite lifetime (via the uncertainty principle) and detector effects.

The relativistic Breit–Wigner distribution is a standard idealized model for this resonance line‑shape.

2.2 Typical real‑world use cases

• Resonance modeling: \(Z^0\), \(W\), \(\phi\), \(\rho\) mesons, …

• Peak finding and mass/width estimation from invariant mass histograms

• Likelihood-based inference in particle physics analyses

2.3 Relations to other distributions

• Cauchy / Breit–Wigner: the classic (non‑relativistic) Breit–Wigner is the Cauchy distribution.

• For large \(\rho\), rel_breitwigner looks locally like a Cauchy near its mode (see Section 5).

• Heavy tails imply that not all moments exist: mean and variance exist, but skewness and kurtosis diverge.

3) Formal Definition

SciPy’s standardized PDF is

\[ f(x\mid \rho) = \frac{k(\rho)}{(x^2 - \rho^2)^2 + \rho^2},\qquad x\ge 0,\ \rho>0. \]

where

\[ k(\rho)=\frac{2\sqrt{2}\,\rho^2\sqrt{\rho^2+1}}{\pi\,\sqrt{\rho^2 + \rho\sqrt{\rho^2+1}}}. \]

CDF

The CDF is the integral

\[F(x\mid\rho)=\int_0^x f(t\mid\rho)\,dt,\]

and has a closed form that SciPy evaluates using complex arithmetic (it’s real‑valued after taking an imaginary part).

Location/scale

With location and scale:

\[ f_Y(y\mid\rho,\text{loc},\text{scale}) = \frac{1}{\text{scale}}\,f\!\left(\frac{y-\text{loc}}{\text{scale}}\,\Big|\,\rho\right), \qquad y\ge \text{loc},\ \text{scale}>0. \]

def _check_positive(name: str, value: float) -> float:
    value = float(value)
    if not np.isfinite(value) or value <= 0:
        raise ValueError(f"{name} must be a finite positive number")
    return value


def rel_breitwigner_k(rho: float) -> float:
    """Normalization constant k(rho) for the standardized PDF."""
    rho = _check_positive("rho", rho)
    return (
        2
        * np.sqrt(2)
        * rho**2
        * np.sqrt(rho**2 + 1)
        / (np.pi * np.sqrt(rho**2 + rho * np.sqrt(rho**2 + 1)))
    )


def rel_breitwigner_pdf(
    x: np.ndarray | float,
    rho: float,
    *,
    loc: float = 0.0,
    scale: float = 1.0,
) -> np.ndarray:
    """PDF of rel_breitwigner with optional loc/scale (NumPy-only)."""
    rho = _check_positive("rho", rho)
    scale = _check_positive("scale", scale)

    x = np.asarray(x, dtype=float)
    y = (x - float(loc)) / scale

    # A numerically stable form used in SciPy: pdf(x) = (k/rho**2) / ( ((x^2-rho^2)/rho)**2 + 1 )
    C = (
        np.sqrt(2 * (1 + 1 / rho**2) / (1 + np.sqrt(1 + 1 / rho**2)))
        * 2
        / np.pi
    )
    denom = (((y - rho) * (y + rho) / rho) ** 2 + 1)
    pdf = (C / denom) / scale

    return np.where(y >= 0, pdf, 0.0)


def rel_breitwigner_logpdf(
    x: np.ndarray | float,
    rho: float,
    *,
    loc: float = 0.0,
    scale: float = 1.0,
) -> np.ndarray:
    """Log-PDF of rel_breitwigner with optional loc/scale (NumPy-only)."""
    rho = _check_positive("rho", rho)
    scale = _check_positive("scale", scale)

    x = np.asarray(x, dtype=float)
    y = (x - float(loc)) / scale

    C = (
        np.sqrt(2 * (1 + 1 / rho**2) / (1 + np.sqrt(1 + 1 / rho**2)))
        * 2
        / np.pi
    )

    denom = (((y - rho) * (y + rho) / rho) ** 2 + 1)
    logpdf = np.log(C) - np.log(denom) - np.log(scale)

    return np.where(y >= 0, logpdf, -np.inf)


def rel_breitwigner_cdf(
    x: np.ndarray | float,
    rho: float,
    *,
    loc: float = 0.0,
    scale: float = 1.0,
) -> np.ndarray:
    """CDF of rel_breitwigner with optional loc/scale (NumPy-only).

    This mirrors SciPy's closed-form implementation, which uses complex arithmetic.
    """
    rho = _check_positive("rho", rho)
    scale = _check_positive("scale", scale)

    x = np.asarray(x, dtype=float)
    y = (x - float(loc)) / scale

    C = np.sqrt(2 / (1 + np.sqrt(1 + 1 / rho**2))) / np.pi

    y_pos = np.where(y >= 0, y, 0.0)
    z = np.sqrt(-1 + 1j / rho) * np.arctan(y_pos / np.sqrt(-rho * (rho + 1j)))
    cdf = C * 2 * np.imag(z)

    cdf = np.clip(cdf, 0.0, 1.0)
    return np.where(y >= 0, cdf, 0.0)


# Quick numerical sanity checks
rho_test = 2.0
xs = np.linspace(0, 50, 200_001)
area = np.trapz(rel_breitwigner_pdf(xs, rho_test), xs)
print("Approx integral of PDF over [0, 50]:", area)
print("CDF(0) =", rel_breitwigner_cdf(0.0, rho_test))
print("CDF(50) =", rel_breitwigner_cdf(50.0, rho_test))

Approx integral of PDF over [0, 50]: 0.9999926082589535
CDF(0) = 0.0
CDF(50) = 0.9999926082589539

4) Moments & Properties

4.1 Mean, variance, skewness, kurtosis

Because the tail behaves like \(f(x)\sim k/x^4\), moments exist only up to order \(n<3\):

\[\mathbb{E}[X^n] < \infty \iff n < 3.\]

For \(X\sim \mathrm{RelBW}(\rho)\) (standardized):

• Mean \(\mathbb{E}[X]\) exists and is finite.

• Variance \(\mathrm{Var}(X)\) exists and is finite.

• Skewness and kurtosis are infinite (SciPy reports nan).

4.2 MGF and characteristic function

• The MGF \(M_X(t)=\mathbb{E}[e^{tX}]\) is finite for \(t<0\) and diverges for \(t>0\) (polynomial tails cannot compete with \(e^{tX}\)).

• The characteristic function \(\varphi_X(t)=\mathbb{E}[e^{itX}]\) always exists.

4.3 Entropy

The differential entropy

\[h(X)=-\int_0^\infty f(x)\log f(x)\,dx\]

does not have a commonly used simple closed form; SciPy evaluates it numerically.

def rel_breitwigner_mean(rho: float) -> float:
    rho = _check_positive("rho", rho)
    A = np.sqrt(rho**2 + 1)
    return (A / np.pi) * np.sqrt(2 * rho / (rho + A)) * (np.pi / 2 + np.arctan(rho))


def rel_breitwigner_second_moment(rho: float) -> float:
    rho = _check_positive("rho", rho)
    return rho * np.sqrt(rho**2 + 1)


def rel_breitwigner_var(rho: float) -> float:
    m = rel_breitwigner_mean(rho)
    return rel_breitwigner_second_moment(rho) - m**2


for rho in [0.2, 1.0, 3.5]:
    m_formula = rel_breitwigner_mean(rho)
    v_formula = rel_breitwigner_var(rho)
    m_scipy, v_scipy = rel_breitwigner.stats(rho, moments="mv")

    print(f"rho={rho}")
    print("  mean  (formula, SciPy):", m_formula, m_scipy)
    print("  var   (formula, SciPy):", v_formula, v_scipy)

print("entropy(rho=1.0) via SciPy:", rel_breitwigner.entropy(1.0))

rho=0.2
  mean  (formula, SciPy): 0.3286859777676976 0.3286859777676975
  var   (formula, SciPy): 0.09592630856260402 0.09592630856260406
rho=1.0
  mean  (formula, SciPy): 0.9653913793583742 0.965391379358374
  var   (formula, SciPy): 0.4822330470336309 0.4822330470336309
rho=3.5
  mean  (formula, SciPy): 3.284899599858133 3.2848995998581327
  var   (formula, SciPy): 1.9496269250927831 1.9496269250927867
entropy(rho=1.0) via SciPy: 0.7677961497524766

def rel_breitwigner_cf(t: float, rho: float) -> complex:
    """Characteristic function via numerical integration (SciPy quad)."""
    rho = _check_positive("rho", rho)
    t = float(t)

    f_cos = lambda x: np.cos(t * x) * rel_breitwigner_pdf(x, rho)
    f_sin = lambda x: np.sin(t * x) * rel_breitwigner_pdf(x, rho)

    re, _ = integrate.quad(f_cos, 0, np.inf, limit=300)
    im, _ = integrate.quad(f_sin, 0, np.inf, limit=300)
    return re + 1j * im


for t in [0.0, 0.5, 1.0, 2.0]:
    val = rel_breitwigner_cf(t, rho=2.0)
    print(f"phi({t}) = {val}")

phi(0.0) = (0.9999999999999997+0j)
phi(0.5) = (0.5630438239678969+0.7174512887274822j)
phi(1.0) = (-0.1597740266239293+0.7104281135014987j)
phi(2.0) = (-0.2863681210052288-0.20255826100750324j)

5) Parameter Interpretation

5.1 Meaning of parameters

In the physics motivation, for a resonance with characteristic mass \(M_0\) and decay width \(\Gamma\) (in natural units), SciPy’s parametrization uses

• shape: \(\rho = M_0/\Gamma\)

• scale: \(\text{scale}=\Gamma\)

• location: typically \(\text{loc}=0\)

so that the mode occurs at

\[\text{mode}(Y)=\text{loc}+\text{scale}\,\rho = M_0.\]

5.2 How the shape changes

• Increasing \(\rho\) moves the peak to the right and makes the distribution more “Cauchy‑like” near its mode.

• Increasing scale stretches the distribution (wider peak; larger absolute variance).

def plot_pdf_cdf_for_rhos(rhos, x_max=10.0):
    xs = np.linspace(0, x_max, 900)

    fig = make_subplots(
        rows=1,
        cols=2,
        subplot_titles=("PDF", "CDF"),
        horizontal_spacing=0.12,
    )

    for rho in rhos:
        fig.add_trace(
            go.Scatter(
                x=xs,
                y=rel_breitwigner_pdf(xs, rho),
                mode="lines",
                name=f"rho={rho}",
            ),
            row=1,
            col=1,
        )
        fig.add_trace(
            go.Scatter(
                x=xs,
                y=rel_breitwigner_cdf(xs, rho),
                mode="lines",
                name=f"rho={rho}",
                showlegend=False,
            ),
            row=1,
            col=2,
        )

    fig.update_xaxes(title_text="x", row=1, col=1)
    fig.update_xaxes(title_text="x", row=1, col=2)
    fig.update_yaxes(title_text="f(x)", row=1, col=1)
    fig.update_yaxes(title_text="F(x)", row=1, col=2)
    fig.update_layout(title="rel_breitwigner: PDF and CDF (standardized)")
    return fig


plot_pdf_cdf_for_rhos([0.5, 1.0, 2.0, 5.0], x_max=12)

# Local Cauchy approximation for large rho near the mode.
# For large rho, around x pprox rho, rel_breitwigner(x; rho) pprox Cauchy(loc=rho, scale=1/2).

rho = 10.0
xs = np.linspace(max(0, rho - 3), rho + 3, 800)

pdf_rbw = rel_breitwigner_pdf(xs, rho)
pdf_cauchy = cauchy.pdf(xs, loc=rho, scale=0.5)

fig = go.Figure()
fig.add_trace(go.Scatter(x=xs, y=pdf_rbw, mode="lines", name="rel_breitwigner"))
fig.add_trace(
    go.Scatter(
        x=xs,
        y=pdf_cauchy,
        mode="lines",
        name="Cauchy approx",
        line=dict(dash="dash"),
    )
)
fig.update_layout(
    title="Large-rho local approximation near the mode",
    xaxis_title="x",
    yaxis_title="density",
)
fig

6) Derivations

6.1 Expectation (mean)

Start from

\[\mathbb{E}[X]=\int_0^\infty x\,\frac{k}{(x^2-\rho^2)^2+\rho^2}\,dx.\]

Substitute \(u=x^2\) so that \(du=2x\,dx\):

\[ \mathbb{E}[X]=\frac{k}{2}\int_0^\infty \frac{1}{(u-\rho^2)^2+\rho^2}\,du. \]

This is a shifted Cauchy kernel. Using

\[\int \frac{1}{(u-a)^2+b^2}\,du = \frac{1}{b}\arctan\left(\frac{u-a}{b}\right),\]

with \(a=\rho^2\) and \(b=\rho\) gives

\[ \mathbb{E}[X] =\frac{k}{2\rho}\left[\arctan\left(\frac{u-\rho^2}{\rho}\right)\right]_{u=0}^{u=\infty} =\frac{k}{2\rho}\left(\frac{\pi}{2}+\arctan(\rho)\right). \]

6.2 Second moment and variance

The second raw moment is

\[\mathbb{E}[X^2]=\int_0^\infty x^2 f(x)\,dx.\]

This integral can be evaluated in closed form (e.g. by extending to an even integral over \(\mathbb{R}\) and applying residue calculus), yielding a remarkably simple result:

\[\mathbb{E}[X^2]=\rho\sqrt{\rho^2+1}.\]

Therefore

\[\mathrm{Var}(X)=\rho\sqrt{\rho^2+1} - \big(\mathbb{E}[X]\big)^2.\]

Higher moments diverge because \(f(x)\sim k/x^4\).

6.3 Likelihood

Given i.i.d. samples \(x_1,\dots,x_n\) from the standardized model (support \(x_i\ge 0\)), the log-likelihood for \(\rho\) is

\[ \log L(\rho\mid x_{1:n}) = n\log k(\rho) - \sum_{i=1}^n \log\left((x_i^2-\rho^2)^2 + \rho^2\right). \]

With location/scale, let \(y_i=(x_i-\text{loc})/\text{scale}\) and add the Jacobian term \(-n\log(\text{scale})\).

def nll_rho_only(rho: float, data: np.ndarray) -> float:
    """Negative log-likelihood for rho (standardized: loc=0, scale=1)."""
    rho = float(rho)
    if rho <= 0:
        return np.inf

    x = np.asarray(data, dtype=float)
    if np.any(x < 0):
        return np.inf

    k = rel_breitwigner_k(rho)
    ll = len(x) * np.log(k) - np.sum(np.log((x**2 - rho**2) ** 2 + rho**2))
    return -ll


# Synthetic MLE demo
rho_true = 2.5
x = rel_breitwigner.rvs(rho_true, size=5_000, random_state=rng)

res = optimize.minimize_scalar(
    nll_rho_only,
    bounds=(1e-3, 50.0),
    args=(x,),
    method="bounded",
)

print("true rho:", rho_true)
print("MLE rho :", res.x)

true rho: 2.5
MLE rho : 2.4911323685576505

7) Sampling & Simulation (NumPy-only)

SciPy can sample from rel_breitwigner, but here we implement a from-scratch sampler that uses only:

• a source of uniforms \(U\sim\mathrm{Uniform}(0,1)\) (NumPy RNG)

• the CDF \(F(x\mid\rho)\)

• numerical inversion (bisection)

Algorithm (inverse CDF via bisection)

For each \(u\in(0,1)\):

1. Find a bracket \([\ell, h]\) with \(F(\ell)\le u \le F(h)\) (start at \(\ell=0\) and expand \(h\) by doubling).

2. Repeat bisection for a fixed number of iterations.

This is robust and easy to implement for any monotone CDF.

def rel_breitwigner_ppf_bisect(
    u: np.ndarray | float,
    rho: float,
    *,
    tol: float = 1e-12,
    max_iter: int = 80,
) -> np.ndarray:
    """Inverse CDF (PPF) via bisection (NumPy-only)."""
    rho = _check_positive("rho", rho)

    u = np.asarray(u, dtype=float)
    if np.any(~np.isfinite(u)):
        raise ValueError("u must be finite")

    out = np.empty_like(u)

    mask0 = u <= 0
    mask1 = u >= 1
    mask = ~(mask0 | mask1)

    out[mask0] = 0.0
    out[mask1] = np.inf

    if not np.any(mask):
        return out

    uu = u[mask]

    lo = np.zeros_like(uu)
    hi = np.full_like(uu, max(1.0, rho))

    # Expand hi until all targets are bracketed
    for _ in range(100):
        cdf_hi = rel_breitwigner_cdf(hi, rho)
        need = cdf_hi < uu
        if not np.any(need):
            break
        hi[need] *= 2
    else:
        raise RuntimeError("Failed to bracket all u values")

    # Bisection
    for _ in range(max_iter):
        mid = 0.5 * (lo + hi)
        cdf_mid = rel_breitwigner_cdf(mid, rho)
        go_right = cdf_mid < uu
        lo[go_right] = mid[go_right]
        hi[~go_right] = mid[~go_right]

        if np.max(hi - lo) < tol:
            break

    out[mask] = 0.5 * (lo + hi)
    return out


def rel_breitwigner_rvs_numpy(rho: float, size: int, *, rng: np.random.Generator) -> np.ndarray:
    u = rng.random(size)
    return rel_breitwigner_ppf_bisect(u, rho)


# Quick check: sample moments vs theory
rho = 2.0
x_samp = rel_breitwigner_rvs_numpy(rho, 30_000, rng=rng)

print("sample mean  :", x_samp.mean())
print("theory mean  :", rel_breitwigner_mean(rho))
print("sample var   :", x_samp.var())
print("theory var   :", rel_breitwigner_var(rho))

sample mean  : 1.8575346488861357
theory mean  : 1.8521891046099694
sample var   : 0.8987076883363847
theory var   : 1.041531475763699

8) Visualization

We’ll visualize:

• the PDF and CDF for different \(\rho\)

• a Monte Carlo histogram compared to the theoretical PDF

rho = 2.0
xs = np.linspace(0, 15, 1200)

samples = rel_breitwigner_rvs_numpy(rho, 50_000, rng=rng)

fig = make_subplots(
    rows=1,
    cols=2,
    subplot_titles=("PDF + Monte Carlo histogram", "CDF"),
    horizontal_spacing=0.12,
)

fig.add_trace(
    go.Histogram(
        x=samples,
        nbinsx=120,
        histnorm="probability density",
        name="samples",
        opacity=0.55,
    ),
    row=1,
    col=1,
)

fig.add_trace(
    go.Scatter(
        x=xs,
        y=rel_breitwigner_pdf(xs, rho),
        mode="lines",
        name="theoretical pdf",
        line=dict(width=2),
    ),
    row=1,
    col=1,
)

fig.add_trace(
    go.Scatter(
        x=xs,
        y=rel_breitwigner_cdf(xs, rho),
        mode="lines",
        name="cdf",
    ),
    row=1,
    col=2,
)

fig.update_xaxes(title_text="x", row=1, col=1)
fig.update_xaxes(title_text="x", row=1, col=2)
fig.update_yaxes(title_text="density", row=1, col=1)
fig.update_yaxes(title_text="probability", row=1, col=2)
fig.update_layout(title=f"rel_breitwigner (rho={rho})")
fig

9) SciPy Integration (scipy.stats.rel_breitwigner)

SciPy exposes the distribution as scipy.stats.rel_breitwigner.

Key methods:

• pdf, logpdf

• cdf, ppf (numerical inversion)

• rvs

• fit

Physical tip: when fitting resonance masses, it is usually appropriate to fix loc=0 (see SciPy docs).

rho = 2.0

x0 = np.array([0.0, 0.5, 1.0, 2.0, 5.0])
print("pdf:", rel_breitwigner.pdf(x0, rho))
print("cdf:", rel_breitwigner.cdf(x0, rho))

s = rel_breitwigner.rvs(rho, size=5, random_state=rng)
print("rvs:", s)

# Fit demo (standardized, so loc should be ~0 and scale ~1)
data = rel_breitwigner.rvs(2.5, size=5_000, random_state=rng)
rho_hat, loc_hat, scale_hat = rel_breitwigner.fit(data, floc=0)
print("fit rho, loc, scale:", rho_hat, loc_hat, scale_hat)

pdf: [0.138329 0.153167 0.212814 0.691646 0.006217]
cdf: [0.       0.071569 0.160358 0.583635 0.99095 ]
rvs: [2.292835 3.588495 0.863896 2.080456 0.086615]

fit rho, loc, scale: 2.5159002307382545 0 1.0010439852719006

# Example: Z0 resonance parameters (Particle Data Group numbers in SciPy docs)
M0 = 91.1876
Gamma = 2.4952
rho = M0 / Gamma

xs = np.linspace(70, 110, 1200)

pdf = rel_breitwigner.pdf(xs, rho, loc=0, scale=Gamma)

fig = go.Figure()
fig.add_trace(go.Scatter(x=xs, y=pdf, mode="lines", name="Z0 model"))
fig.update_layout(
    title="Relativistic Breit–Wigner example (Z0 line shape)",
    xaxis_title="invariant mass M (GeV)",
    yaxis_title="density",
)
fig

10) Statistical Use Cases

10.1 Hypothesis testing (signal vs background)

A common pattern is testing for the presence of a resonant signal on top of a smoother background.

We’ll do a toy likelihood‑ratio test between:

• \(H_0\): background only (uniform)

• \(H_1\): mixture of background + rel_breitwigner signal (known parameters)

10.2 Bayesian modeling

You can treat rel_breitwigner as a likelihood for \(\rho\) (and possibly scale), assign priors, and compute a posterior.

10.3 Generative modeling

In simulation pipelines, rel_breitwigner acts as a parametric generator of resonance masses.

def logpdf_uniform_bg(x: np.ndarray, L: float) -> np.ndarray:
    x = np.asarray(x, dtype=float)
    return np.where((0 <= x) & (x <= L), -np.log(L), -np.inf)


def logpdf_signal(x: np.ndarray, rho: float, scale: float) -> np.ndarray:
    return rel_breitwigner_logpdf(x, rho, loc=0.0, scale=scale)


def logpdf_mixture(x: np.ndarray, w: float, rho: float, scale: float, L: float) -> np.ndarray:
    # log( w*s + (1-w)*b ) computed stably
    x = np.asarray(x, dtype=float)
    log_s = logpdf_signal(x, rho, scale)
    log_b = logpdf_uniform_bg(x, L)

    a = np.log(w) + log_s
    b = np.log(1 - w) + log_b

    m = np.maximum(a, b)
    return m + np.log(np.exp(a - m) + np.exp(b - m))


# Toy LRT
rng = np.random.default_rng(7)

L = 200.0
rho = 36.5
scale = 2.5
w = 0.15
n = 400

# observed data under H1
x_sig = rel_breitwigner.rvs(rho, loc=0, scale=scale, size=int(n * w), random_state=rng)
x_bg = rng.uniform(0, L, size=n - len(x_sig))
x_obs = np.concatenate([x_sig, x_bg])

ll0 = np.sum(logpdf_uniform_bg(x_obs, L))
ll1 = np.sum(logpdf_mixture(x_obs, w=w, rho=rho, scale=scale, L=L))
T_obs = 2 * (ll1 - ll0)
print("Observed LRT statistic T:", T_obs)

# Monte Carlo null distribution under H0
B = 400
T_null = np.empty(B)
for b in range(B):
    x0 = rng.uniform(0, L, size=n)
    ll0_b = np.sum(logpdf_uniform_bg(x0, L))
    ll1_b = np.sum(logpdf_mixture(x0, w=w, rho=rho, scale=scale, L=L))
    T_null[b] = 2 * (ll1_b - ll0_b)

p_value = np.mean(T_null >= T_obs)
print("Monte Carlo p-value:", p_value)

fig = go.Figure()
fig.add_trace(go.Histogram(x=T_null, nbinsx=40, name="T under H0"))
fig.add_vline(x=T_obs, line_width=2, line_dash="dash", line_color="black")
fig.update_layout(
    title="Toy likelihood-ratio test: background vs (background+signal)",
    xaxis_title="T",
    yaxis_title="count",
)
fig

Observed LRT statistic T: 96.15435313231501
Monte Carlo p-value: 0.0

# Simple Bayesian inference for rho (standardized: scale=1, loc=0)
# Posterior p(rho | x) ∝ p(x | rho) p(rho)

rng = np.random.default_rng(7)

rho_true = 2.0
x = rel_breitwigner.rvs(rho_true, size=400, random_state=rng)

rho_grid = np.linspace(0.2, 8.0, 600)

# Prior: log-normal on rho (weakly informative)
mu, sigma = 0.0, 0.7
log_prior = -np.log(rho_grid * sigma * np.sqrt(2 * np.pi)) - (np.log(rho_grid) - mu) ** 2 / (2 * sigma**2)

log_like = np.array([-nll_rho_only(r, x) for r in rho_grid])
log_post_unnorm = log_like + log_prior
log_post = log_post_unnorm - np.max(log_post_unnorm)
post = np.exp(log_post)
post = post / np.trapz(post, rho_grid)

rho_map = rho_grid[np.argmax(post)]

print("true rho:", rho_true)
print("MAP rho :", rho_map)

fig = go.Figure()
fig.add_trace(go.Scatter(x=rho_grid, y=post, mode="lines", name="posterior"))
fig.add_vline(x=rho_true, line_width=2, line_dash="dash", line_color="black")
fig.add_vline(x=rho_map, line_width=2, line_dash="dot", line_color="gray")
fig.update_layout(
    title="Posterior over rho (scale fixed to 1)",
    xaxis_title="rho",
    yaxis_title="density",
)
fig

true rho: 2.0
MAP rho : 2.0230383973288815

# Generative modeling: draw synthetic resonance masses (Z0 example)

rng = np.random.default_rng(7)

M0 = 91.1876
Gamma = 2.4952
rho = M0 / Gamma

masses = Gamma * rel_breitwigner_rvs_numpy(rho, 30_000, rng=rng)

fig = go.Figure()
fig.add_trace(
    go.Histogram(
        x=masses,
        nbinsx=140,
        name="simulated events",
        histnorm="probability density",
        opacity=0.7,
    )
)

xs = np.linspace(70, 110, 900)
fig.add_trace(
    go.Scatter(
        x=xs,
        y=rel_breitwigner.pdf(xs, rho, loc=0, scale=Gamma),
        mode="lines",
        name="model pdf",
    )
)

fig.update_layout(
    title="Generative simulation: Z0 resonance masses",
    xaxis_title="invariant mass M (GeV)",
    yaxis_title="density",
)
fig

11) Pitfalls

• Invalid parameters: rho and scale must be strictly positive.

• Support matters: in standardized form, \(x<0\) has density 0; with location, support is \(x\ge\text{loc}\).

• Higher moments diverge: sample skewness/kurtosis can be unstable or meaningless.

• Fitting loc: allowing loc to vary can produce non-physical negative mass support; in resonance modeling, prefer floc=0.

• Numerical CDF/PPF: SciPy’s ppf relies on numerical inversion; for extreme probabilities it can be slow. For simulation, rvs is usually preferable.

12) Summary

• rel_breitwigner is a continuous distribution on \([0,\infty)\) used to model relativistic resonance line shapes.

• PDF: \(f(x\mid\rho)=\dfrac{k(\rho)}{(x^2-\rho^2)^2+\rho^2}\) with \(\rho>0\).

• Mean and variance exist, but skewness/kurtosis diverge.

• Sampling can be done via numerical inverse CDF (NumPy-only) or via SciPy’s rvs.

• In physics parametrization: \(\rho=M_0/\Gamma\) and scale=Γ, so the mode is at \(M_0\).